import numpy as np
import datetime

import volsurface
import math
import normdist

from mscf.fc3.project.volfitter.analytics.volsurface import  ParametricVolSurface

STRIKE_BY_SPOT_LOWER_BOUND = 0.5
STRIKE_BY_SPOT_UPPER_BOUND = 1.5
NUMBER_OF_STRIKES = 100

NUM_DAYS_IN_YEAR = 365.

class VolSurfaceFitter :

    def __init__(self, optionMktData):
        self.optionMD=optionMktData

    def fitParametricVolSurface(self) :
        chainMapKeys=self.optionMD.optionPriceChainMap.keys()
        chainMapValues=self.optionMD.optionPriceChainMap.values()

        ##      pDate=self.optionMD.pricingDate
        keyCount=len(chainMapKeys)

        mMatrix=[]

        for index in range(len(chainMapKeys)):
            params=self.fitParametricCurve(chainMapKeys[index],chainMapValues[index])
            mMatrix.append(params)
        ##            print len(params)
        ##            for i in range(len(params)):
        ##                print params[i]


        expDates = sorted(self.optionMD.optionPriceChainMap.keys())
        today = datetime.date.today()
        terms = [((expDate - self.optionMD.pricingDate).days / NUM_DAYS_IN_YEAR)  for expDate in expDates]

        dtms=[]
        for expDate in expDates:
            term = (expDate - today).days / NUM_DAYS_IN_YEAR
            dtm = volsurface.DTM(expDate, term)
            dtms.append(dtm)

        return ParametricVolSurface(self.optionMD.refSpot, mMatrix, self.optionMD.interestRate, terms, dtms)

    def fitParametricCurve(self,  expDate, priceChain ):
    ##        print "fitParametricCurve for one priceChain--Same Term--different Strike"
        T=((expDate-priceChain.pricingDate).days) /NUM_DAYS_IN_YEAR
        ##        pricePairMap=priceChain.getCallPutPairMap()
        pricePairs=priceChain.getCallPutPairMap().values()
        callPriceArray = priceChain.callMktPriceArray
        ##        print pricePairs

        for index in range(len(callPriceArray)):
#            midPrice = pricePairs[index].callPrice.midPrice
#            if not midPrice is None:
#                pricePairs[index].callPrice.impVol=self.ImpliedVolatility(self.optionMD.refSpot,pricePairs[index].strike,self.optionMD.interestRate, T,pricePairs[index].callPrice.midPrice)
            callPrice = callPriceArray[index]
            midPrice = callPrice.midPrice
            if not midPrice is None:
                callPrice.impVol=self.ImpliedVolatility(self.optionMD.refSpot,callPrice.strike,self.optionMD.interestRate, T,callPrice.midPrice)

        ##        for index in range(len(pricePairs)):
        ##            print pricePairs[index].callPrice.midPrice
        ##            print pricePairs[index].callPrice.impVol

        L=len(callPriceArray)
        S=self.optionMD.refSpot
        r=self.optionMD.interestRate

        F = S * math.exp(r * T)
        Vatm = 1

        for i in range(L-1):
            if callPriceArray[i].strike < F and callPriceArray[i+1].strike >= F:
                callPriceArray
                Vatm = (callPriceArray[i+1].impVol - callPriceArray[i].impVol)/(callPriceArray[i+1].strike-callPriceArray[i].strike) * (F - callPriceArray[i].strike) + callPriceArray[i].impVol
                #Vatm = (pricePairs[i+1].callPrice.impVol-pricePairs[i].callPrice.impVol)/(pricePairs[i+1].strike-pricePairs[i].strike) * (F -pricePairs[i].strike) + pricePairs[i].callPrice.impVol
                break

        mMatrix = [[0 for i in range(3)] for j in range(L)]
        Vol = [[0 for i in range(1)] for j in range(L)]

        for i in range(L):
            mMatrix[i][0] = 1
            mMatrix[i][1] = math.log(callPriceArray[i].strike/F)/(Vatm*math.sqrt(T))
            mMatrix[i][2] = mMatrix[i][1] * mMatrix[i][1]
            Vol[i][0] = callPriceArray[i].impVol

        A = np.matrix(mMatrix)
        V = np.matrix(Vol)
        B = A.I

        parameters = B*V
        qFit = [[0 for x in range(2)] for y in range(L)]

        for i in range(L):
            qFit[i][0] = callPriceArray[i].strike
            qFit[i][1] = parameters[0] + mMatrix[i][1]*parameters[1] + mMatrix[i][2]*parameters[2]
        ##                print qFit[i][0]
        ##                print pricePairs[i].callPrice.impVol
        ##                print qFit[i][1]

        paramsWithVatm=[T,S,r,Vatm, parameters[0],parameters[1],parameters[2] ]
        return paramsWithVatm

    @staticmethod
    def ImpliedVolatility( S, X, RF, T,MPR,MAXERR=0.08, optionType="call", YIELD=0 ):
        #S = stock price
        #X = exercise price
        #RF = risk-free rate
        #T = time to expiration
        #MPR = market price of option
        #MAXERR = maximum acceptable difference between market price and model price
        #YIELD = dividend yield
        #optionType--Call, put

        Pi=3.1415926535;

        sigma=math.sqrt(abs(math.log(S * math.exp(-YIELD * T) / X) + RF * T) * (2 / T))
        impVol=sigma

        Diff=10000000

        while Diff > MAXERR * MPR:

            d1=(math.log(S * math.exp(-YIELD * T) / X) + (RF + sigma*sigma / 2) * T) / (sigma * math.sqrt(T))
            d2 = d1 - sigma * math.sqrt(T)

            Nd1=normdist.phi(d1)
            Nd2=normdist.phi(d2)
            if optionType=="call":
                BSP = S * math.exp(-YIELD * T) * Nd1 - X * math.exp(-RF * T) * Nd2
            else:
                BSP = X * math.exp(-RF * T) * (1 - Nd2) - S * math.exp(-YIELD * T) * (1 - Nd1)

            impVol=sigma
            sigma = sigma - ((BSP - MPR) * math.exp(-d1 * d1 / 2) * math.sqrt(2 * Pi)) / (S * math.exp(-YIELD * T) * math.sqrt(T))

            Diff=abs(MPR-BSP)

        return impVol;

##
##        spot=108.99
##        interestRate=0.0047        



##        iMaturityCount=len(optionMarketData.optionPriceChainMap)
##        strikeLowerBound = spot * STRIKE_BY_SPOT_LOWER_BOUND
##        strikeUpperBound = spot * STRIKE_BY_SPOT_UPPER_BOUND
##
##        mMatrix = [[0 for i in range(3)] for j in range(NUMBER_OF_STRIKES)]
##
##
##        #parameters = B*V
##
##        paramVolSurface = ParametricVolSurface(mMatrix, parameters)
##
##        return paramVolSurface
